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ABSTRACT 

Many astrophysical models predict a diffuse flux of high-energy neutrinos from active galactic nuclei 
and other extra-galactic sources. At muon energies above 1 TeV, the upward-going muon flux induced 
by neutrinos from active galactic nuclei is expected to exceed the flux due to atmospheric neutrinos. 
We have performed a search for this astrophysical neutrino flux by looking for upward-going muons 
in the highest energy data sample from the Super-Kamiokande detector using 1679.6 live days of 
data. We found one extremely high energy upward-going muon event, compared with an expected 
atmospheric neutrino background of 0.46 ±0.23 events. Using this result, we set an upper limit on the 
diffuse flux of upward-going muons due to neutrinos from astrophysical sources in the muon energy 
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range 3.16-100 TeV. 

Subject headings: galaxies: active — gamma rays: bursts — neutrinos 



1. INTRODUCTION 

The GeV-PeV energy range is unexplored territory for 
neutrino astronomy — observations of neutrinos at these 
energies will open a new window on the high energy uni- 
verse. A wide variety of astrophysical phenomena are 
expected to produce extremely high energy neutrinos, 
ranging from activ e galactic nuclei (AGNs) and gamma- 
ray b ursts (GRBs: iHalzen fc HoopeHr2002t iGaisser et all 
|1995() to more exotic sources such as dark matter anni- 
hilation or decays of topological defects (Stasto 2004i). 

The flux of neutrinos at such high energies is quite 
small; therefore, large-scale detectors are required. One 
effective technique for observing high-energy neutrinos 
with an underground detector is to look for muons 
produced by i/^ or interacting in the surrounding 
rock. (Throughout this paper, "muons" will refer to 
both ^+ and The muon range in rock increases 

with muon energy, which expands the effective interac- 
tion volume for high-energy events. Downward-going 
neutrino-induced muons cannot be distinguished from 
the much larger flux of downward cosmic-ray muons, but 
since cosmic ray muons cannot travel through the entire 
Earth, upward-going muons are almost always neutrino- 
induced. Thus, upward-going muons provide a suitable 
high-energy neutrino signal. 

At muon energies above 1 — 10 TeV, the upward- 
going muon flux due to neutrinos from AGNs is ex- 
pect ed to exceed the flux du e to atmospheric neutri- 
nos (jStecker fc SalamonI Il996t iMannheim et al.l |2001D . 
This cosmic neutrino flux could be detected either by 
searching for point sources of high-energy neutrinos or 
by detecting a diffuse, isotropic flux of neutrinos from 
unresolved astrophysical sources. A diffuse cosmic neu- 
trino flux would be observed as an excess to the ex- 
pected atmospheric neutrino flux at high energies. In 
this analysis, we focus on searching for a diffuse flux of 
upward-going muons due to neutrinos from astrophysical 
sources using the highest energy data sample in Super- 
Kamiokande (Super-K). This study complements other 
Supcr-K searches for astrophysical point sources of high 
energy neu trinos that use data over a larger energy range 
()Abe et al.i 2006). In this paper we describe a search for 
evidence of a high energy astrophysical neutrino flux in 
Super-K's highest energy upward-going muon sample. In 
§ [2] we describe the Super-Kamiokande detector, and in 
§[3] we give the details of how we selected candidate events 
from Super-K's ultra-high-energy sample. We evaluate 
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our selection process with Monte Carlo in § [4] and calcu- 
late the observed upward-going muon flux in § [5l Sec- 
tions [H] and [7] discuss the background due to the atmo- 
spheric neutrino flux. Based on the results, we set an 
upper limit in § [8] and conclude in § [9l Any necessary 
estimates and approximations have been made so that 
they lead to a conservative result for this upper limit. 

2. THE SUPER-KAMIOKANDE DETECTOR 

The Super-K detector is a cylindrical 50 kiloton wa- 
ter Cerenkov detector, located in the Kamioka-Mozumi 
mine in Japan. It is 41.4 m tall and 39.3 m in diameter. 
The detector was constructed under the peak of Mount 
Ikenoyama, which provides an average rock overburden 
of 1000 m (2700 m water equivalent). Its geodetic loca- 
tion is at 36.4° north, 137.3° east, and altitude 370 m. 

Super-K consists of two concentric, optically separated 
detectors. Until 2001 July the inner detector (ID) was 
instrumented with 11,146 inward- facing 50 cm diameter 
photomultiplier tubes (PMTs). The outer detector (OD) 
is a cylindrical shell of water surrounding the ID and is 
instrumented with 1885 outward-facing 20 cm diameter 
PMTs. Between the ID and the OD, there is a 50 cm 
thick shell. Photons coming from this region will not be 
detected by either the OD or the ID, so we refer to it as 
the insensitive region. 

More detai l s abo ut the detector can be found in 
iFukuda et al.l (|2003D . The data sample used in this anal- 
ysis was taken from 1996 April to 2001 July, correspond- 
ing to 1679.6 days of detector livetime. This data run is 
referred to as SK-I. 

Super-K is primarily designed to detect lower energy 
neutrinos from the Sun, the atmosphere, and particle ac- 
celerators but can potentially detect the extremely high 
energy neutrinos expected from astrophysical sources as 
well. This paper focuses on the events at the highest 
energy end of Super-K's detection range. 

3. EVENT SELECTION 

The ultra-high-energy sample in SK-I consists of 
events that deposit > 1.75 x 10^ photoelectrons (pe) in 
the ID. In the low-energy regime, on average about 9 pe 
are recorded by the ID PMTs for each MeV of energy de- 
posited in the tank; the electronics for the ID PMTs sat- 
urate at about 300 pe. Thus an event with > 1.75 x 10^ 
pe in the ID corresponds to a minimum energy deposition 
of approximately 200 GeV, but the actual energy depo- 
sition could be much higher, since the saturation effect 
prevents all of the produced pe from being recorded. 

At high energies, muons have some probability 
to lose energy through radiative processes such as 
bremsstrahlung, resulting in an electromagnetic shower 
that deposits large quantities of pe in the detector. For 
comparison, a muon that traverses the maximum path 
length through the ID (50 m) but does not produce any 
electromagnetic showers will deposit approximately 11 
GeV via ionization energy loss, corresponding to ^ 10^ 
pe deposited in the ID. Thus a high-pe cutoff offers a 
means of selecting high energy events. 

At the high-pe threshold of > 1.75 x 10'^ pe, the 
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high level of saturation in the ID PMT electronics can 
cause Super-K's precision muon fitting algorithms to fail. 
Therefore, these extremely energetic events are not in- 
cluded in other studies of upward-going muons in SK- I 
(|Desai et alJl2004l : iFukuda et aOl999l : lAshie et al.l[2005f ) . 
In this study we analyzed this ultra-high-energy data 
sample separately using a different fitting method based 
on information from the OD. 

3.1. Outer Detector Linear Fit 

SK-I's ultra-high-energy data sample contains a to- 
tal of 52214 events. Most of these are either very en- 
ergetic downward-going cosmic-ray muons or multiple 
muon events where two or more downward-going muons 
hit the detector simultaneously. In order to select candi- 
date upward-going muons from this sample, we applied 
a simple linear fit to the OD data for each event. A 
linear fit was done on the z-position of each OD PMT 
versus the time it fired, weighted by the total charge in 
the PMT. Example fits of simulated downward-going and 
upward-going muon events are shown in Figure [1] The 
slope of this fitted line is an estimate of — cos 8, where 
is the zenith angle of the incoming muon. A positive 
fitted slope indicates that the muon is upward-going. A 
similar linear fit was done on the x- and y-positions to 
determine the full muon trajectory through the detector. 

Since this fitting method is based on the OD (which 
has a lower resolution than the ID), it is not as precise 
as the muon fitting algorithms used in the lower energy 
upward-going muon analysis. However, it works even 
when the ID PMT electronics are completely saturated 
and the precision ID-based algorithms fail. 

3.2. Selection Cuts 

To select candidate upward-going muons, we applied 
the OD-based fit to all 52214 events in the ultra-high- 
energy data sample. A cut of cos 6 < 0.1 was used to 
eliminate the bulk of the downward-going single and mul- 
tiple muon events. In addition, the fitted trajectory was 
required to have a path length of > 7 m in the ID. 

To ensure high-quality fit results, we looked at the 
number of OD PMTs hit near the projected entry and 
exit points. For a true throughgoing muon, there will 
be a cluster of hit PMTs around the entry and exit 
points. If the fit is accurate, the projected path should 
pass through both of these clusters, so we made an addi- 
tional cut on the number of OD PMTs hit within 10 m 
of the projected OD entry and exit points, iVoDentry and 
A^ODoxit- We required both A^oDcntry and A^oDoxit to be 
10 or greater. 

Events that do not have A^oDcntry and A^oDexit > 10 
are generally either stopping muons, partially contained 
events, or poorly fitted throughgoing muons. Stopping 
muons are muons that stop in the detector and only form 
an OD entry cluster. They typically have energies of 1-10 
GeV, well below the range of the expected astrophysical 
signal, so it is appropriate to discard events that look 
like stopping muons. Partially-contained events are neu- 
trino interactions that take place inside the detector and 
only form an OD exit cluster — they are not part of the 
upward-going muon fiux incident on the detector, so we 
want to discard these as well. This cut also occasionally 
eliminates inaccurately fitted throughgoing muon events. 



which reduces the efficiency somewhat but improves the 
accuracy of the fit results. 

Another possible type of event that can masquerade as 
a throughgoing muon is a partially-contained event with 
multiple exiting particles. Such an event will create two 
(or more) clusters of hit PMTs in the OD, which could be 
mistaken as the entry and exit points of a throughgoing 
muon. In order to eliminate such events, we looked at 
the timing between the OD and ID entry points. If the 
event is truly a throughgoing muon, the OD PMTs near 
the entry point should fire before the ID PMTs. We 
determined OD and ID entry clusters via a simple time- 
based clustering method and evaluated the mean time of 
hits within 6 m of the OD and ID entry points, tioentry 

and toDontry 

In an ideal measurement, tiDcntry < ^ODontry would in- 
dicate that the perceived entry cluster is actually caused 
by an exiting particle. However, the timing determina- 
tion is complicated by an effect known as prepulsing. 
Prepulsing occurs when there are so many photons in- 
cident on the PMT that they are not all converted to 
photoelectrons at the photocathode - some photons hit 
the first dynode instead and are converted to photoelec- 
trons there. When this happens in an ID PMT, the PMT 
pulse appears earlier relative to the OD light. If this oc- 
curs for an ultra- high energy throughgoing muon event, 
it could cause tiDentry " ^ODentry to be negative. To al- 
low for this effect, we used a fairly loose cut of —40 ns: 
if tiDontry " toDontry < "40 US, indicating very early ID 
light, then the event was rejected as a likely neutrino 
interaction in the ID. 

After these cuts on cos 9, path length, AoDentry, 
A^ODcxit, and iiocntry -toDcntry wcrc appficd to the 52214 
events in the sample, 343 candidate events remained. 
These remaining events were then evaluated by a visual 
scan and a manual direction fit by two independent re- 
searchers to select events with cos 8 < 0. 

The visual scan eliminates events that can pass the au- 
tomatic reduction but are obvious to the trained human 
eye as noise, i.e., mainly downward-going multiple muon 
events and "fiashers." Multiple muon background events 
occur when two or more downward-going muons from an 
atmospheric shower pass through the detector simultane- 
ously. They have extra energy deposition (and thus are 
expected to be more common in the ultra-high-energy 
sample) and typically give poor OD fit results due to 
multiple OD clusters, but they are easy to identify visu- 
ally. Flashers are events caused by malfunctioning PMTs 
that emit light and create characteristic patterns. After 
the multiple muon events and flashers are removed, the 
manual direction fit separates truly upward-going muons 
from mis-fitted downward-going muons. 

From the 343 candidates, only one event passed the 
visual scan and manual direction fit selection as being 
truly upward-going. The breakdown of the visual scan 
and manual fit classifications is shown in Tabic [T] 

This upward-going muon event selected from the > 
1.75 X 10^ pe sample is the ultra- high energy upward- 
going muon signal observed by SK-I. This event occurred 
on 2000 May 12 at 12:28:07 UT and deposited 1,804,716 
pe in the ID. Based on the manual fit results, the path 
length through the ID was 40 m, and the zenith angle was 
cos© = —0.63, corresponding to a direction of origin of 
(R.A., deck) = 20^38™, -37°18'. 



4 



Swanson et al. 




TABLE 1 

Visual scan of candidate upward-going muon events 



Visual scan classification 


Number of events 


Multiple muon events 


164 


PMT "flashers" (malfunctioning PMTs) 


5 


Other noise events 


2 


Downgoing muons (manual fit cos > 0) 
Upgoing muons (manual fit cos © < 0) 


171 
1 


Total 


343 



4. HIGH-ENERGY ISOTROPIC MONTE CARLO 

In order to calculate the observed muon flux, we need 
to determine the resolution and efhciency of the OD- 
based fit and other cuts on high-energy muons, and we 
need to estimate the probability that a muon of a given 
energy will deposit > 1.75 x 10^ pe in the ID. 

To determine these quantities, we generated a high- 
energy isotropic Monte Carlo (MC) sample. This MC 
consists of an isotropic flux of muons in monoenergetic 
bins impinging on the Super-K detector, representing a 
flux of muons from neutrino interactions in the surround- 
ing rock. Seven monoenergetic bins were used, with 
muon energies ranging from 100 GeV to 100 TeV, and 
10000 events with a path length in the ID of > 7 m 
were generated in each bin. The simulation was per- 
formed using a GEANT-based detector simulation. GEANT's 
muon propagation has been shown to agree with the- 
oretical pr edictions up to muon energies of 100 TeV 
(jBottai fc Pcrrone 2001; Desai et al. 2003). 

4.1. Resolution and Efficiency of Event Selection 

The muon trajectory fitting algorithm discussed in 
§ 13.11 was applied to the high-energy isotropic MC. A 
plot illustrating the angular resolution of the fit is shown 
in Figure [21 The resolution of 5-12° is poor compared 
to the typical 1° resolution of precision ID fitting algo- 
rithms, but those algorithms do not work well on these 
high-energy saturated events. 

The efficiency of the upward-going cuts was estimated 
by considering all of the MC events with true values of 
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angle between fitted and true direction [degrees] 

Fig. 2. — Angular resolution of OD-based fit for events from the 
high energy isotropic MC with > 10^ pe in the ID. 

cos O < and ID path length > 7 m and then determin- 
ing the fraction of these events that pass the selection 
cuts described in § 13.21 This was done using isotropic 
MC events with > 10^ pe since there were relatively few 
events in the MC with > 1.75 x 10^ pe. (Above 10^ pe the 
efficiency and resolution do not depend strongly on the 
number of ID pe deposited.) Also, as discussed in § 14. 2( 
only the energy bins in the range 3.16 — 100 TeV were 
considered. The efficiency was calculated as a function 
of muon energy and cos Q . 

Statistical uncertainties on the efficiency determina- 
tio n were ca l culate d using the Bayesian method discussed 
bv lConwavl (120021 ). Systematic uncertainties due to un- 
certainties in the fitted values of the zenith angle, path 
length, and number of OD PMTs hit near the entry and 
exit points were evaluated by comparing the fitted val- 
ues with the true MC values. The fitted values agree well 
with the MC values, and distributions of the difference 
between the fitted values and the MC values were used 
to determine an effective 1 a region for each parameter. 
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The cuts on these parameters were varied by 1 cr in either 
direction to determine the effect on the efficiency. 

Another source of systematic uncertainty on the effi- 
ciency is prepulsing, which is not included in the MC sim- 
ulation and must be estimated separately. To do this, we 
compared the results from the high-energy isotropic MC 
to a sample of 627 ultra-high-energy downward-going 
muon data events with > 1.75 x 10^ pe in the ID se- 
lected by a visual scan. For many of these data events, 
the time difference tiocntry — ioDcntry is negative, showing 
evidence of prepulsing not seen in the MC sample. 

To estimate the systematic uncertainty on the effi- 
ciency calculation due to prepulsing, we calculated a 
heuristic correction to the MC by smearing the MC dis- 
tribution of tiDontry — ^ODontiy such that it matched the 
distribution of the ultra-high-energy downward-going 
data. We applied this adjustment to the events in the 
MC sample and found that an additional 1% of the 
MC events would be cut after accounting for prepuls- 
ing, which lengthens the negative-side error bars on the 
efficiency estimates by approximately 0.01. 

Finally, one additional correction must be made: the 
efficiency has been estimated at different values for iJ^, 
but the flux calculation is for the flux above a thresh- 
old energy i?™™. Since the efficiency decreases with 
energy, e(i?^,6) is an overestimate for e(>i?™'",8). 
This does not lead to a conservative upper limit for the 
flux, so a correction must be made. This requires knowl- 
edge of the energy spectrum of the expected signal, so 
we modeled the signal as an isotropic flux of neutri- 
nos with d^i, (E^)/dE^ oc E ~'^ (a plausible astrophys- 
ical spectrum iGaissed Il990l l , used the method of § [7] 
to estimate the muon flux 'I'p (> -B™'"), and used this 
muon flux to extrapolate our efficiency calculations to 
find e (> i?™'",9). This procedure yields small down- 
ward adjustments (0.6 — 9%) to the calculated efficiency 
in each angular bin. 

Additionally, we must also account for the efficiency of 
the visual scan and manual fit procedure. To do this, we 
did a visual scan of the 605 events from the high-energy 
isotropic MC with true cos < 0, true ID path length 
> 7 m, and > 1.75 x 10^ pe in the ID and found that 
10 of these events were eliminated in the visual scan. 
This gives an efficiency of roughly 98%, so we adjust our 
efficiency results by a factor of 0.98. 

The final results for the efficiency of our cuts on up- 
ward, throughgoing, > 1.75 x 10® pe muons (including 
all corrections discussed above) are plotted in Figure [3l 

4.2. Ultra-High- Energy Fraction 

In order to set an upper limit on the flux of upward- 
going muons from cosmic neutrinos one must make an 
inference about the energies of the upward-going muons 
in the ultra-high-energy sample. Above energies of 
~ 1 TeV, muon energy loss in water is dominated by ra- 
diative processes such as bremsstrahlung, so high-energy 
muons have some probability of depositing large numbers 
of photoelectrons in the Super-K detector and contribut- 
ing to the ultra-high-energy sample. Since this energy 
loss is not continuous, it is not possible to estimate the 
muon energy for a single ultra-high-energy upward-going 
muon event. Rather, MC is used to make a statistical 
statement about the energies of the muons that make up 
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Fig. 3. — Efficiency of our data reduction on upward, through- 
going: ^ 1.75 X 10® pe muons as determined by the high-energy 
isotropic MC. Error bars include statistical and systematic errors. 



TABLE 2 

Fraction of high-pe events k(E^) in high energy isotropic 

MC 



number of MC events k {Ei_t) 





(out of 10000) with 


(with statistical 


(TeV) 


> 1.75 x 10® pe in ID 


uncertainties) 


0.1, 0.316, 1 





0.0000 

-0.0000 


3.16 


35 


0.0035 ± 0.0006 


10 


113 


0113 +°-°°^^ 
^■^^^^'^ -0.0010 


31.6 


293 


0.0293 ± 0.0017 


100 


879 





the > 1.75 X 10® pe sample. 

The high-energy isotropic MC has been used to deter- 
mine the fraction k {E^j) of muons with energy E^ that 
will deposit > 1.75 x 10® pe in the ID, thus contributing 
to the ultra-high-energy sample. Results are shown in 
Table [21 Statistical uncertainties were calculated using 
the Bayesian method discussed bv IConwavl (|2002l ). As 
can be seen in Tabled the three lowest energy bins — 
from 100 GeV to 1 TeV — do not make a significant 
contribution to the > 1.75 x 10® pe sample. Hence, the 
rest of this analysis was done using only the four highest 
energy bins — from 3.16 to 100 TeV. 

These results for k [E^) are used in §[5]to calculate the 
flux $p (> -E™'"). Since k{Efj) increases with energy, 
it is an underestimate for /c (> £^™™), which leads to a 
conservative upper limit for (> £^™"). 

Also, since this MC only simulates the muon and not 
the actual neutrino interaction, the effect of lower energy 
debris from deep inelastic scattering events that make it 
into the detector has been neglected, again leading to 
an underestimate of k{E^). This effect is expected to 
be small in the energy range considered here — above 
Ef^ = 1 TeV at the detector entry, over 80% of the 
neutrino-induced upward-going muons come from over 
200 m away from the detector, so most of the debris is 
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TABLE 3 

The flux of ultra-high energy upward-going muons as 
observed by sk-i. 



(TeV) (cm-'s-lsr-l) 



3.16 


2.64 X 10" 


-14 


-H6.1% 
-17.9% 


10 


8.23 X 10" 


-15 


-(-9.48% 
-9.73% 


31.6 


3.25 X 10~ 


-15 


-1-6.67% 
-6.40% 


100 


1.10 X 10" 


-15 


-1-4.96% 
-4.09% 



Note. — These fluxes include both atmospheric background 
and potential astrophysical signal at each threshold energy. 

absorbed by the surrounding rock, as determined using 
the atmospheric neutrino MC discussed in § [6l 

5. FLUX CALCULATION 

The flux of upward-going muons above a threshold en- 
ergy is given by 

n . 

(j) ^min\ _ \ ^ 

mV- p ) - 2T,Tk{>Ef^)j^^e{>Ef^,Q,)A{Q, 

where n is the total number of upward-going muon events 
observed and Qj is the zenith angle of the jth event. The 
efficiency £ (> -E™'", Qj) and the ultra-high-energy frac- 
tion k (> -E™'") are calculated in § 21 T is the detector 
livetime, which is 1679.6 days for SK-I. A{Qj) is the 
effective area of the Super-K detector perpendicular to 
the direction of incidence for tracks with a path length 
of > 7 m in the ID. The average effective area of the 
detector is ^ 1200 m^. 

Equation ([T]) has been applied to the detected upward- 
going muon event discussed in § 13.21 to calculate 
*M (> -E™™) for in the range 3.16 - 100 TeV. Re- 

sults are shown in Table [3l Systematic uncertainties 
include a 0.1% uncertainty on the live time T, a 0.3% 
uncertainty on the effective area A, the total efficiency 
uncertainties shown in Figure [31 and the statistical un- 
certainties on k shown in TablejH This flux includes both 
the potential signal from astrophysical neutrinos and a 
background of atmospheric neutrinos. 

6. EXPECTED ATMOSPHERIC BACKGROUND 
FROM MONTE CARLO 

When searching for neutrinos from astrophysical 
sources, the dominant background is the atmospheric 
neutrino spectrum. Atmospheric neutrinos are produced 
by decays of pions and kaons formed when cosmic rays 
interact with particles in the atmosphere. We have used 
an atmospheric neutrino MC that is a 100 yr equivalent 
sample of events du e to the atmospher ic neutrino flux. 
The neutrino flux in lHonda et al.l (|2004( ) was used up to 
ne utrino energies of 1 TeV. At 1 TeV, the calculated flux 
m iVolkoval (fT98l was rescaled to the Honda et al. flux. 
Above 1 TeV, the rescaled flux from Volkova was used up 
to 100 TeV. Neutrino interactions wer e modeled using th e 
GRV94 parton distribution functions (jGluck et alJIlQQSl ) , 
and muon propagation through the rock and water was 



modeled using GEANT. Furthe r details on the atm ospheric 
neutrino MC can be found in lAshie et al.l (|2005[ ). No cor- 
rection is made for neutrino oscillations, be cause based 
on the oscillation parameters determined in lAshie et al.l 
()2005f ). the neutrino oscillation probability is negligible 
for neutrinos above 1 TeV. 

This atmospheric MC is split into two parts: 
a partially-contained/fully-contained (PC/FC) sample, 
which consists of events with neutrino interaction points 
inside the ID plus a shell 50 cm thick surrounding the 
ID (the insensitive region), and an upward-going muon 
sample, which consists of events with neutrino interac- 
tion points outside the ID. Note that these two samples 
overlap because they both cover the 50 cm insensitive 
region. 

The OD-based fit was applied to the events in the at- 
mospheric MC, using the same cuts that were applied 
to the SK-I data. A total of 11 MC events passed 
the > 1.75 X 10*5 pe, cos 9 < 0.1, path length > 7 m, 
iVoDcntry and -/VoDoxit ^ 10, OD/ID timing, and manual 
fit cuts. Out of these 11, 2 are from the PC/FC sample, 
both with interaction points inside the ID. The remain- 
ing 9 events are from the upward-going muon sample: 3 
events with interaction points in the 50 cm insensitive 
region, 1 event in the water of the OD, and 5 events in 
the rock surrounding the detector. 

All of these background events are deep inelastic scat- 
tering (DIS) events where an interaction between a muon 
neutrino and a nucleon produces a muon plus a spray of 
lower energy particles. The 6 events with interaction 
points within the detector (ID or OD) have muon ener- 
gies of 0.1 — 0.8 TeV, and the 5 events occurring in the 
rock have muon energies of 1 — 20 TeV. This difference in 
the energy range can be understood as follows: For DIS 
events occurring a long distance (> 2 m or so) from the 
detector, only the muon will reach the detector since the 
lower energy debris will be absorbed by the surrounding 
rock, but for nearby events or events occurring in the 
water of the OD, some of these lower energy particles 
will enter the detector as well. This means that nearby 
events can be included in the > 1.75 x 10^ pe sample 
with lower muon energies than more distant events. 

Since the insensitive region is covered by both the 
PC/FC and the upward-going muon MC samples, we 
divided the 3 events originating from this region in half, 
for a total of 1.5 events in the insensitive region. This 
gives a total of 9.5 MC events in 100 yr of simulated live 
time. Scahng the 100 yr MC to SK-I's live time of 1679.6 
days gives an expected background of 0.44 events due to 
atmospheric neutrinos during the operation of SK-I. 

The statistical uncertainty in this background mea- 
surement of the MC events is 31%. There are also signif- 
icant systematic uncertainties: the normalization of the 
atmospheric neutrino flux has a theoretical uncertainty 
of ±1 0% at neutrino energies below 10 GeV (|Ashie et al.l 
l2005f ). In order to extend this to the energy range of 
the expected background, we must also account for the 
uncertainty of 0.05 in the spectral index of the primary 
cosmic ray spectrum above 100 GeV, which leads to a 
0.05 uncertainty in the spectral index for atmospheric 
neutrinos above 10 GeV ()Ashie et al.ll2005f ). 

To determine how much this uncertainty affects our 
result for the background, we consider the atmospheric 
neutrino flux to be known at 10 GeV, and we calculate 
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TABLE 4 

Systematic uncertainties in atmospheric neutrino 
background 



Source of Uncertainty 


Uncertainty 


Statistical 


31% 


Absolute normalization of atm v flux 


10% 


Primary spectral index 


37% 


Neutrino cross-section uncertainty 


10% 


Total uncertainty in background flux 


50% 



amount of energy in the detector. In contrast, the peak 
90% of the expected muon events from a harder spectrum 
— a hypothetical astrophysical flux — come from 
the range E^, = l - 6000 TeV. Thus, even though the 
atmospheric flux is very small in the energy range = 
3 — 100 TeV and above where we are setting our limit, 
the high-pe tails of the distribution from lower energy 
events dominate our background simply because of the 
much larger flux of lower energy atmospheric events. 



the uncertainty of the total flux above a threshold 
energy 10.6 TeV, the average neutrino energy of the at- 
mospheric MC events passing our cuts. For a differential 
flux of d^„/dE^ (X E^'^ with 7 = 3.7 ±0.05, the spectral 
index uncertainty gives us a ±37% uncertainty on the 
atmospheric neutrino flux 

Finally, the neutrino cross-section at high energies is 
thought to be known to within 10% or less, so we include 
an additional 10% uncertainty to account for this. These 
uncertainties are summarized in Table [?] and lead to a 
total uncertainty on the background of 50%. 

Other potential errors (uncertainties in the simulation 
of the SK detector and varying hadron multiplicities in 
different deep inelastic scattering models) were tested 
and shown to not make a significant contribution to the 
systematic uncertainty. 

Another potential background source of high-energy 
neutrinos not included in the 100 yr atmospheric MC is 
the prompt atmospheric neutrino flux, which arises from 
decays of short-lived charmed particles produced when 
cosmic rays interact with particles in the atmosphere. 
This flux is not as well-understood as the conventional at- 
mospheric neutrino flux from decays of pions and kaons, 
but it is expected to have a harder spectrum and there- 
fore is expected to become more important as we push 
towards higher energy scales. 

Based on the 100 yr atmospheric MC, we calculate 
an expected background for this analysis of 0.44 ± 0.22 
events, compared to the 1 event observed. However, 
there are three effects that this MC does not take into ac- 
count: it does not include neutrinos over 100 TeV, it does 
not include the prompt atmospheric neutrino flux, and 
it does not account for attenuation of neutrinos passing 
through the Earth. We account for these issues by mak- 
ing corrections based on an analytical calculation dis- 
cussed in § [T] 

Also, it is important to note that the atmospheric back- 
ground — both conventional and prompt — comes from 
a lower energy range than that for which we expect to 
observe the possible signal of neutrinos from astrophys- 
ical sources. Roughly speaking, the peak 90% of the 
expected upward-going muon events come from muons 
with energies in the range Ef^ = 0.02 — 10 TeV for con- 
ventional atmospheric neutrinos, and E^^ = 0.2—200 TeV 
for prompt neutrinos. As we learned from the 11 back- 
ground events in the 100 yr atmospheric MC, muons with 
energies Ef^ < 1 TeV contribute to the 1.75 x 10^ pe sam- 
ple mainly via debris from DIS events very close to the 
detector rather than from catastrophic energy loss of the 
muon. 

Since there are many more low-energy events in the 
atmospheric spectrum, they will dominate even though 
each one only has a tiny probability of depositing a large 



7. ANALYTICAL ESTIMATE OF EXPECTED 
MUON FLUX 

7.1. Method for Calculating Muon Flux 

In order to better understand our observed flux of high- 
energy upward-going muons, we developed a method to 
calculate the expected upward-going muon event rate due 
to a predicted flux of neutrinos. We have used this to plot 
curves for theoretical muon fluxes in Figure [5] and also 
to adjust the atmospheric background calculated with 
MC in § [6] by correcting for effects not included in the 
simulation. 

To convert a model neutrino flux into an expected 
upward-going muon event rate, we follo w the calcula- 
tion detailed in iGaisser et al.l ()1995f ) and I Gandhi et al.l 
|l996|). The flux of muons 

*M (> -^m'") above an energy 

threshold ii"™'" is given by 



*M (> E, 



min \ 
M ) 



dE^P, (E^ 



rpmin 



dE^ ' 



where Pu(E^, 



(2) 



E'r^ 



is the probability that an incom- 
ing neutrino with energy E^, will produce a muon with 
energy above the threshold -E™™ at the detector, and 
c?<I>^^ {E^)/dEi, is the differential neutrino flux averaged 
over solid angle and reduced by an exponential factor due 
to attenuation of the neutrinos as they pass through the 
Earth. 

Pfj, {Eu, -E™™) depends on both the charged-current 
neutrino cross-section and the energy lost by the re- 
sulting muon as it propagates through the Earth. In 
our calculation, we calculate the cross-se ction using the 
GRV94 parton distribution fu nctions ( Cluck et al.l 19950 
apphed with the code used in lGandhi et al.l (|1996 ) pro- 
vided by M. Reno (2005, private communication), and we 
use an effective muon rang e calculated using MC meth- 
ods (|Lipari fc Stane v "1991V 

As discussed in ^andhi e t"al] ()1996f ). the appropriate 
cross-section a {E,^) to use in the neutrino attenuation 
factor lies between the charged current cross section and 
the sum of the charged and neutral current cross sections, 
so we have calculated the flux with the upper and lower 
limits and averaged the resulting fluxes together. 

We integrated equation ^ numerically using various 
predicted models of the neutrino flux d^^/dEi, to calcu- 
late theoretical (> i;™™) curves. These are shown in 
Figure [5] for comparison to our experimental limits. 

7.2. Analytical Estimate of Expected Background 

Of particular interest here are the model fluxes of the 
background due to atmospheric neutrinos, both conven- 
tional and prompt. We use the method discussed in ii l7.1l 
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■ 100 year aim MC 
- HE isotropic MC 
Power law fit to HE isotropic MC 



10"' 



ll I 1^1 III I 

10^ 10^ 10" 10^ 

Muon Energy [GeV] 

Fig. 4. — Fraction of events k{E^) that deposit > 1.75 X 10® pe 
in the ID for the 100 yr atmospheric MC {open circles), the high- 
energy isotropic MC {filled squares), and a power-law fit to the re- 
sults from high-energy isotropic MC. The power law is a reasonable 
approximation for the dominant energy range of the atmospheric 
flux {E^ = 100 GeV - 30 TeV). 

to correct for the omissions in the 100 yr atmospheric 
MC discussed in § |6] by accounting for neutrinos over 
100 TeV, attenuation of neutrinos in the Earth, and the 
flux of prompt neutrinos. 

The expected number of events N seen by Super-K in 
livetime T is given by 



TABLE 5 

Confidence intervals for the upward-going muon flux due 
to neutrinos from astrophysical sources. 





90% C.L. ran 


ge 


(TeV) 


( cm~^s~-'-sr~ 


') 


3.16 


- 1.03 X 10" 


13 


10 


- 3.19 X 10" 


14 


31.6 


- 1.26 X 10- 


14 


100 


0-4.28 X 10" 


15 



We also used this calculation to correct for the prompt 
atmospheric neutrino flux. To account for the theoretical 
uncertainty in the prompt flux due to differences between 
various flux models, we defined a high model and a low 
model for the prompt flux that bra cket the models shown 
in Figure 1 of lGelmini et al.l (j2003[ ) that are not ruled out 
by experimental lim its. The models used a re discussed 
in mo r e detail in iThurmian et al.l (I1996D. IZas et al.l 
11993'), 'Rvazhskava et al.l (I2002D. iBu gaev et al.' ('1989'). 
IPasquah et al. (1999I). and rCelmini et al. (2000a,b). Our 
analytical calculation gives -/Vioo = 0.033 events for the 
low model and iVioo = 0.94 events for the high model, 
corresponding to a prompt flux that is 6.7% ± 6.2% of 
the flux of conventional atmospheric neutrinos. 

Based on these results, we correct the background esti- 
mate made using the 100 yr atmospheric MC in §[6] by ap- 
plying these relative scalings to the result of 0.44 events 
in the SK-I livetime from § [6l This gives us a final re- 
sult for the background of 0.46 ±0.23 events for the SK-I 
exposure. 



TV = 2TrTA, 



d<i>^ (> Ef-) 



(3) 



)/dE, 



is the derivative of the curve 
"""^ is the 



5'"'" that will de- 
We estimated k ( ET^'^ 



where d^f, 

calculated by the method in § 17.1! and k [E\ 
fraction of muons with energy above E^^ 

posit > 1.75 X 10^ pe in the ID. 

using the results of the high-energy isotropic MC dis- 
cussed in ii l4.2l bv fitting the data in Table[2]with a simple 
power law with a cutoff at 1 since k represents a prob- 
ability. Figure S] compares this power-law fit to results 
from 100 yr atmospheric MC, illustrating that the power 
law is reasonable in the energy range we are considering. 

We analytically estimated the expected number of 
events in the 100 yr atmospheric MC by starting with 
the same input neutrino fiux, applying equation ^ with- 
out the exponential neutrino attenuation factor and in- 
tegrating up to a maximum neutrino energy of E^^'^^ = 
100 TeV. We then used this muon flux in equation ^ 
to find the expected number of background events in 100 
yr and obtained A^ioo — 7.39 expected events. (The 100 
subscript denotes 100 yr of exposure.) This matches our 
results from the 100 yr atmospheric MC to within sta- 
tistical uncertainties. 

To estimate the effects of neutrinos above 100 TeV, 
we repeated the above calcuation without the E^'^^ = 
100 TeV cutoff and obtained A^ioo = 7.47 events, an in- 
crease of 1.1%. Including the neutrino attenuation factor 
as well, we obtained A^ioo = 7.29 ±0.02 events per 100 yr 
of exposure, a decrease of 2.4% ±0.3%. (The uncertainty 
is due to the choice of cross-section.) 



8. UPPER LIMIT FOR MUON FLUX FROM 
COSMIC NEUTRINOS 

Using the observed ultra-high-energy upward-going 
muon signal of 1 event and the expected atmospheric 
neutrino background of 0.46 ± 0.23 events, we have cal- 
culated 90% confidence upper limits for the upward-going 
muon fiux in the 3.16 — 100 TeV range due to neutrinos 
from astrophysical sources (or any other non-atmospheric 
sources). 

Thi s was done using the method of lFeldman fc Cousins! 
(|1998l ). with the sy stem atic un certaintie s incor porated 
using the met hod of ICou sins & Highland! (|1992t l. as im- 
plemented bv ' Conrad et al.i (j200§) and improved bv lHilll 
|2003). This method incorporates both uncertainties in 
the background fiux and uncertainties in the flux factor / 
relating the observed number of events n to the observed 
flux <I>: ^ ~ f n. The uncertainty in / includes system- 
atic errors in the livetime, effective area, efficiency, and 
ultra-high-encrgy fraction. For the confldence interval 
calculation, the largest percent error for each energy bin 
from Table [3] was used as the percent uncertainty in /. 
The uncertainties in both the background and the flux 
factor were assumed to have a Gaussian distribution. 

The final results are shown in Table [5] and plot- 
ted in Figure [5l along wit h models of variou s pos- 
sible signals from A GNs (jStecker fc SalamonI 119961; 
Man nheim et al.ll2001l ) and GRBs (jWaxman fc Bahcalll 
1 19991) . as well as the backgrounds due to atmospheric 
neutrinos (Honda et al^2004; Volkova 1980) and prompt 
neutrinos (Gclmini et al.ll2003i ). 

The upper limits calculated here are consistent with 
the models of astrophysical signals. Also shown are 
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TABLE 6 

Approximate upper limits from SK-I on astrophysical 

NEUTRINOS {Vfj, + Ufj). 



I Atmospheric flux 

I Prompt flux 
AGN core model (SS) 
AGN jet model (IVIPR) 
GRB model (WB) 



[GeV] 



Fig. 5. — Upper limits from this analysis on muon (/i''' + /i ) 
flux above energy threshold E™'", compared to various model 
fluxes. Models shown for muon flux due to astrophysical neu- 
trinos are AGN models from SS l lStecker fc Salamon 1996) and 
MPR (Mannheim et al.' '20011), and a GRB model from WB 
fWaxman & Bahcall 1999). Also shown is t he atmospheric back- 
ground, as modeled by .Honda et al.l II2004I ) below = I TeV 
and by Volkova (1980]) rescaled to match the Honda model above 
El, = 1 TeV. The upper edge of the atmospheric band represents 
the horizontal flux, and the lower edge represents the vertical flux. 
The background due to muons from prompt atmospheric neutrinos 
(assumed to be is otropic) is shown for a range of possible models 
as summarized in 'Ge lmini et al.l ||2003l) . Finally, we also show the 
expected muon flux from a model neutrino flux that is isotropic 
and proportional to i?,7^ at the values of the limit set by MACRO 
ijAmbrosio et al. 2003) and AMANDA-II I GroB et'al .: 2005;). The 
models of the neutrino flux have been converted to a muon flux with 
equ ation (|2]l using the GRV94 parton distribut ions (Gluck et ahj 
I1995il and the effective muon range from lLipari fc Stanev (1991). 



models with a hypothetical isotropic neutrino flux with 
a spectrum proportional to E^^^ and a normalization 
scaled to the limits o n an flux set b y MACRO 

(lAmbrosio et al.l [200l and AMANDA-II (jCroB et all 
|2005() . The model neutrino fluxes were converted into 
muon fluxes using equation ([2]) as discussed in § [T] 

To facilitate easier comparison with other experiments, 
we also convert our limits on the muon flux into approx- 
imate limits on the neutrino flux. In order to do this, we 
assume a model neutrino flux that is isotropic and pro- 
portional to E~'^ . To get an approximate neutrino limit, 
we find normalization factors for an -E^^ muon flux curve 
such that the curve passes through each of our four limit 
points in Figure O and we use these factors to find the 
implied limits on E^^ flux. 

In order to determine the approximate neutrino energy 
range in which these limits arc valid, we use equation ^ 
to determine the neutrino energy range that produces 
the bulk of the muon signal for a given value of the muon 
energy threshold E™^'^. We define the energy range as the 
range that (1) produces 90% of the muon flux above E™™ 
and (2) has a higher value of the integrand of equation ^ 
within the range than anywhere outside the range. This 
definition is based on the definit ion of hi g hest p osterior 
density intervals as described bv lConwavl (120021 ). 

The results of these approximations for neutrino limits 
and energy ranges are shown in Table [H] and are also 
plotted in Figure [6l with the same models and experi- 



^min 

(TeV) 


90% C.L. upper limit 
(GeV cm^^s^-'-sr^-'-) 


Neutrino energy range 
(GcV) 


3.16 


6.0 X 10-5 


6.3 X 10^ - 1.4 X 10^ 


10 


3.7 X 10^5 


1.7 X 10-* - 2.4 X 10^ 


31.6 


3.2 X 10-5 


4.7 X 10'' - 5.1 X 10** 


100 


2.6 X 10-5 


1.4 X 105 - 1.1 X lO'^ 



Note. — Upper limits on E'^ {d^i, /dE,^). Note that converting 
from a muon flux limit to a neutrino flux limit requires additional 
assumptions — our limits on the muon flux are shown in Table [5] 



SK-I limits 



Atmospheric flux 
Prompt flux 

AGN core model (SS) 

■ AGN jet modef (f^/PR) 
GRB model (WB) 




10° 
Ev [GeV] 



10'" 



Fig. 6. — Approximate upper limits from SK-I on astrophysical 
neutrinos {Vfj, -|- P^). Models shown are the same as in Figure [5] 
Note that converting from a muon flux limit to a neutrino flux limit 
requires additional assumptions - our limits on the muon flux are 
shown in Figure [5] 

mental limits as shown in Figure [S] To draw our limits 
on this plot, we chose to draw only the most sensitive 
limit in each energy range. 

Note that these results for the limits on the neutrino 
flux are only approximations made to facilitate compar- 
ison to other experiments — our primary results arc the 
limits on the muon flux shown in Table and Figure [S] 

9. CONCLUSIONS 

In conclusion, we have developed a method for analyz- 
ing Super-K's highest energy data to search for evidence 
of high-energy neutrino flux from astrophysical sources. 
We have done a thorough study of the efficiency and the 
expected backgrounds from this method and applied our 
method to the SK-I data sample. Our study of the high- 
est energy events in SK-I does not show evidence of a 
high-energy cosmic neutrino signal. 

We have set upper limits on the muon flux due to 
cosmic neutrino sources. These limits arc consistent 
with the results of oth er experiments (Ambrosio et al] 
I2OO3I: iCrofi et all 120051 ). It is possible that an astro- 
physical neutrino signal could be within the grasp of the 
next generation of neutrino detector s such as I ceCube 
()Ahrens et al.l[200l and ANTARES ()Katzll2004( l. 
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